##Supplementary File-1
## Stable Isotope Mixing Models in R with simmr.
## Calculateing the contribution of various moisture sources in Zagros region karstic groundwater samples
## Instal "install.packages('simmr')" at first
  
library(simmr)

##karstic samples in Orumieh station
mix = matrix(c(-9.65, -10.48, -7.79, -10.3, -10.3, -9.72, -10.01, -10.14, -9.55, -9.46, -8.39, -10.41, -10.36, -63.9, -66.2, -60.0, -67.1, -64.8, -64.1, -65.6, -64.8, -63.2, -62.8, -60.0, -66.8, -66.8, 13.3, 17.64, 2.32, 15.3, 17.6, 13.66, 14.48, 16.32, 13.2, 12.88, 7.12, 16.48, 16.08), ncol = 3, nrow = 13)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7, -20.8, -25.4, -19.2, -24.4, 7.3, 16.5, 18.5, 21.6), ncol = 3, nrow = 4)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 25.42, 25.01, 24.38, 27.81, 8.60, 6.19, 6.40, 10.30), ncol = 3, nrow = 4)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)


 simmr_out = simmr_mcmc(simmr_in)
 summary(simmr_out, type = 'diagnostics')
 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')
 
 end 
 
 library(simmr)

##karstic samples in Keyno station

mix = matrix(c(-6.91, -5.83, -5.23, -5.69, -6.16, -4.96, -5.42, -6.14, -5.7, -33.79, -25.16, -25.32, -24.8, -27.22, -20.05, -25.42, -26.47, -23.49, 21.49, 21.48, 16.52, 20.72, 22.06, 19.63, 17.94, 22.65, 22.11), ncol = 3, nrow = 9)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7, -20.8, -25.4, -19.2, -24.4, 7.3, 16.5, 18.5, 21.6), ncol = 3, nrow = 4)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 25.42, 25.01, 24.38, 27.81, 8.60, 6.19, 6.40, 10.30), ncol = 3, nrow = 4)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)


 simmr_out = simmr_mcmc(simmr_in)
 
 
 summary(simmr_out, type = 'diagnostics')
 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')
 
end

library(simmr)

##karstic samples in kangan station

mix = matrix(c(-3.4, -4.35, -14.6, -19.7, 12.6, -15.1), ncol = 3, nrow = 2)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass", "mT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7,-3.8, -20.8, -25.4, -19.2, -24.4,-22.9,7.3, 16.5, 18.5, 21.6, 8.3), ncol = 3, nrow = 5)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 2.2, 25.42, 25.01, 24.38, 27.81, 24.3, 8.60, 6.19, 6.40, 10.30, 9.3), ncol = 3, nrow = 5)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)


 
 simmr_out = simmr_mcmc(simmr_in)
 summary(simmr_out, type = 'diagnostics')
 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')
 
end

library(simmr)

##karstic samples in Hormozgan station

mix = matrix(c(-4.3, -4.76, -5.1, -4.57, -6.5, -4.21, -4.57, -4.47, -5.11, -5.5, -5.36, -22.8, -24.9, -24.2, -23.7, -35.5, -20.7, -22.1, -22.8, -25.8, -29.5, -25.5, 11.6, 13.18, 16.6, 12.86, 16.5, 12.98, 14.46, 12.96, 15.08, 14.5, 17.38), ncol = 3, nrow = 11)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass", "mT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7,-3.8, -20.8, -25.4, -19.2, -24.4,-22.9,7.3, 16.5, 18.5, 21.6, 8.3), ncol = 3, nrow = 5)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 2.2, 25.42, 25.01, 24.38, 27.81, 24.3, 8.60, 6.19, 6.40, 10.30, 9.3), ncol = 3, nrow = 5)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)


 simmr_out = simmr_mcmc(simmr_in)
 summary(simmr_out, type = 'diagnostics')
 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')
 
end

library(simmr)

##karstic samples in Tang Bijar station

mix = matrix(c(-4.12, -4.76, -6.29, -3.85, -3.86, -17.3, -19.5, -28.3, -16.1, -16.0, 15.66, 18.58, 22.02, 14.7, 14.88), ncol = 3, nrow = 5)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7, -20.8, -25.4, -19.2, -24.4, 7.3, 16.5, 18.5, 21.6), ncol = 3, nrow = 4)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 25.42, 25.01, 24.38, 27.81, 8.60, 6.19, 6.40, 10.30), ncol = 3, nrow = 4)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)


 simmr_out = simmr_mcmc(simmr_in)
 summary(simmr_out, type = 'diagnostics')
 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')
 
end


library(simmr)

##karstic samples in Lorestan-Aleshtar

mix = matrix(c(-7.66, -7.61, -7.16, -7.58, -6.85, -6.97, -6.58, -7.75,-40.2, -42.0, -37.1, -41.2, -36.2, -36.4, -33.3, -38.3, 21.08, 18.88, 20.18, 19.44, 18.6, 19.36, 19.34, 23.7 ), ncol = 3, nrow = 8)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7, -20.8, -25.4, -19.2, -24.4, 7.3, 16.5, 18.5, 21.6), ncol = 3, nrow = 4)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 25.42, 25.01, 24.38, 27.81, 8.60, 6.19, 6.40, 10.30), ncol = 3, nrow = 4)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)


 simmr_out = simmr_mcmc(simmr_in)
 summary(simmr_out, type = 'diagnostics')
 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')

 end

library(simmr)

##karstic samples in Kermanshah-Alvand

mix = matrix(c(-6.78, -7.31, -6.45, -5.97, -6.18, -6.63, -6.22, -6.04, -5.95, -35.7, -38.3, -34.8, -32.3, -32.0, -35.0, -34.0, -31.8, -30.7, 18.54, 20.18, 16.8, 15.46, 17.44, 18.04, 15.76, 16.52, 16.9), ncol = 3, nrow = 9)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7, -20.8, -25.4, -19.2, -24.4, 7.3, 16.5, 18.5, 21.6), ncol = 3, nrow = 4)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 25.42, 25.01, 24.38, 27.81, 8.60, 6.19, 6.40, 10.30), ncol = 3, nrow = 4)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)


 simmr_out = simmr_mcmc(simmr_in)
 summary(simmr_out, type = 'diagnostics')
 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')
end

library(simmr)

##karstic samples in Arjan plain

mix = matrix(c(-5.67, -5.7, -5.91, -5.95, -3.52, -4.5, -4.52, -4.82, -2.37, -2.86, -2.63, -5.02, -2.96, -4.25, -3.13, -2.04, -5.05, -4.37, -5.09, -4.68, -4.63, -4.7, -26.4, -25.9, -29.1, -29.8, -13.9, -19.7, -16.4, -22.0, -12.3, -9.1, -9.5, -21.2, -4.8, -17.6, -16.6, -6.0, -20.0, -14.4, -21.4, -17.9, -18.3, -17.2, 18.96, 19.7, 18.18, 17.8, 14.26, 16.3, 19.76, 16.56, 6.66, 13.78, 11.54, 18.96, 18.88, 16.4, 8.94, 10.32, 20.4, 20.56, 19.32, 19.54, 18.74, 20.4), ncol = 3, nrow = 22)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7, -20.8, -25.4, -19.2, -24.4, 7.3, 16.5, 18.5, 21.6), ncol = 3, nrow = 4)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 25.42, 25.01, 24.38, 27.81, 8.60, 6.19, 6.40, 10.30), ncol = 3, nrow = 4)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)

 simmr_out = simmr_mcmc(simmr_in)
 summary(simmr_out, type = 'diagnostics')
 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')
 
end

library(simmr)

##karstic samples in Khersan dam site

mix = matrix(c(-7.22, -6.79, -6.79, -6.84, -6.56, -6.23, -5.9, -6.31, -5.82, -4.92, -6.26, -6.01, -6.83,-6.76, -6.74, -6.73, -6.67, -5.61, -5.87, -5.42, -5.0, -5.62, -6.52, -6.36,-36.95, -35.45, -33.1,  -34.7, -34.65, -32.7, -28.8, -30.05, -27.8, -25.05, -30.1,-29.15, -34.6, -33.252, -33.8, -33.2, -34.5, -30.15, -29.9, -25.3, -24.75, -31.6, -34.05, -31.95, 2.81, 18.87, 21.26, 20.02, 17.87, 17.14, 18.4, 20.47, 18.76, 14.35, 19.98, 18.93, 20.04, 20.83, 20.16, 20.64, 18.86, 14.73, 17.06, 18.1, 15.25, 13.4, 18.15, 18.93), ncol = 3, nrow = 24)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7, -20.8, -25.4, -19.2, -24.4, 7.3, 16.5, 18.5, 21.6), ncol = 3, nrow = 4)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 25.42, 25.01, 24.38, 27.81, 8.60, 6.19, 6.40, 10.30), ncol = 3, nrow = 4)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)


 simmr_out = simmr_mcmc(simmr_in)
 summary(simmr_out, type = 'diagnostics')

 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')
 end

library(simmr)
##karstic samples in Paveh station
mix = matrix(c(-7.30, -7.51, -7.25, -6.55, -7.43, -6.50, -6.78, -7.15, -7.13, -6.91, -6.28, -7.34, -40.9, -40.50, -38.66, -35.50, -36.82, -35.81, -35.89, -37.65, -36.96, -36.31, -33.57, -39.05, 18.28, 19.57, 19.31, 16.90, 22.62, 16.18, 18.36, 19.57, 20.10, 18.96, 16.70, 19.63 ), ncol = 3, nrow = 12)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7, -20.8, -25.4, -19.2, -24.4, 7.3, 16.5, 18.5, 21.6), ncol = 3, nrow = 4)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 25.42, 25.01, 24.38, 27.81, 8.60, 6.19, 6.40, 10.30), ncol = 3, nrow = 4)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)

 simmr_out = simmr_mcmc(simmr_in)
 summary(simmr_out, type = 'diagnostics')

 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')
 end

library(simmr)
##karstic samples in Dareh shahr
mix = matrix(c(-5.87, -6.0, -5.3, -5.6, -5.7, -6.0, -5.7, -6.0, -5.7, -6.8, -6.5, -6.8, -5.8, -5.7, -5.4, -6.2, -5.6, -6.0, -5.6, -5.5, -4.5, -35, -32, -31, -32, -28, -34, -28, -33, -32, -34, -33, -32, -24, -34, -27, -28, -33, -27, -25, 11.4, 16.0, 11.4, 12.8, 17.6, 14.0, 17.6, 21.4, 20.0, 20.4, 13.4, 13.6, 19.2, 15.6, 17.8, 20.0, 11.8, 17.0, 11.0), ncol = 3, nrow = 19)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7, -20.8, -25.4, -19.2, -24.4, 7.3, 16.5, 18.5, 21.6), ncol = 3, nrow = 4)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 25.42, 25.01, 24.38, 27.81, 8.60, 6.19, 6.40, 10.30), ncol = 3, nrow = 4)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)

 simmr_out = simmr_mcmc(simmr_in)
 summary(simmr_out, type = 'diagnostics')

 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')
 end

library(simmr)
##karstic samples in Bakhtegan
mix = matrix(c(-6.76, -6.41, -5.79, -5.78, -6.1, -6.24, -6.85, -6.0, -6.63, -5.87, -7.12, -6.53, -6.22, -6.64, -6.29, -6.7, -4.37, -6.39, -6.82, -6.56, -37.8, -35.1, -33.5, -32.9, -30.8, -32.1, -36.5, -33.2, -33, -35.1, -40.1, -35.7, -33.0, -36.8, -36.1, -37.0, -23.90, -33.1, -37.2, -35.0, -29.0, -27.8, -36.2, 16.28, 16.18, 12.82, 13.34, 18.0, 17.82, 18.3, 14.8, 20.04, 11.14, 16.86, 16.54, 16.76, 16.32, 14.22, 16.6, 11.06, 18.02, 17.36, 17.48, 21.0, 21.32, 22.36, 16.28), ncol = 3, nrow = 24)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7, -20.8, -25.4, -19.2, -24.4, 7.3, 16.5, 18.5, 21.6), ncol = 3, nrow = 4)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 25.42, 25.01, 24.38, 27.81, 8.60, 6.19, 6.40, 10.30), ncol = 3, nrow = 4)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)

 simmr_out = simmr_mcmc(simmr_in)
 summary(simmr_out, type = 'diagnostics')

 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')
 end

library(simmr)
##karstic samples in Turkey
mix = matrix(c(-12.1, -14.0, -12.6, -11.3, -11.5, -10.8, -12.5, -12.8, -12.8, -12.8, -12.9, -11.4, -10.1, -11.3, -11.4, -12.2, -13.5, -11.6, -12.1, -12.7, -10.7, -10.9, -10.7, -9.6, -9.8, -10.0, -87, -96, -94, -91, -92, -79, -95, -90, -94, -91, -83, -79, -82, -90, -85, -91, -85, -83, -89, -82, -72, -69, -65, -66, -69, 9.8, 16.0, 6.8, -0.6, 0, 7.4, 5.0, 12.4, 8.4, 12.2, 8.2, 1.8, 8.4, 1.2, 12.6, 17.0, 7.8, 13.8, 12.6, 3.6, 15.2, 16.6, 11.8, 12.4, 11.0), ncol = 3, nrow = 25)
colnames(mix) = c('d18O','d2H','d-excess')
s_names = c("cP-airmass", "cT-airmass", "mP-airmass", "MedT-airmass")
s_means = matrix(c(-3.5, -5.2, -6.6, -5.7, -20.8, -25.4, -19.2, -24.4, 7.3, 16.5, 18.5, 21.6), ncol = 3, nrow = 4)
s_sds = matrix(c(3.0, 1.74, 4.51, 1.93, 25.42, 25.01, 24.38, 27.81, 8.60, 6.19, 6.40, 10.30), ncol = 3, nrow = 4)
simmr_in = simmr_load(mixtures = mix,
                     source_names = s_names,
                     source_means = s_means,
                     source_sds = s_sds)

 simmr_out = simmr_mcmc(simmr_in)
 summary(simmr_out, type = 'diagnostics')

 summary(simmr_out, type = 'statistics')
 summary(simmr_out, type = 'quantiles')
 end
